%Verlocity-Verlet 算法
% ini_state为一个cell类型变量；t0的单位为us
function state = v_verlet(N, ini_state, t0, steps, K, alpha)
    tau = t0 * 1e-6;
    
    state = ini_state;
    past_state = ini_state;
    
    for i = 1:steps
        for j = 1:N
            r_past = past_state{j}(1);
            v_past = past_state{j}(2);
            f_past = f(N, past_state, j, K, alpha);
            
            r_now = r_past + v_past * tau + 0.5 * f_past * tau^2;
            
            state{j}(1) = r_now;   
        end
        
        for k = 1:N
            v_past = past_state{k}(2);
            
            f_past = f(N, past_state, k, K, alpha);
            f_now = f(N, state, k, K, alpha);
            
            v_now = v_past + 0.5 * (f_past + f_now) * tau;
            
            state{k}(2) = v_now;
        end
        
        past_state = state;
        
    end
    
    return
    
end

%根据构象state计算力
function force = f(N, state, num, K, alpha)

    if(num == 1)
        force = K * (-2 * state{1}(1) + state{2}(1)) - alpha * (state{1}(1)^2 - (state{2}(1) - state{1}(1))^2);
    elseif(num == N)
        force = K * (-2 * state{N}(1) + state{N-1}(1)) - alpha * (-state{N}(1)^2 + (state{N}(1) - state{N-1}(1))^2);
    else
        force = K * (-2 * state{num}(1) + state{num+1}(1) + state{num-1}(1)) - alpha * ((state{num}(1) - state{num-1}(1))^2 - (state{num}(1) - state{num+1}(1))^2);
    end
    
    return
end